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Abstract 

A finite-difference time domain solution of the airfoil gust problem is obtained 
using a high-accuracy nonlinear computational aeroacoustics code. For computational 
efficiency, the equations are cast in chain-rule curvilinear form, and a structured 
multiblock solver is used in parallel. In order to fully investigate the performance of this 
solver, a test matrix of eight problems are computed (two airfoil geometries, two gust 
frequencies, and two gust configurations). These results are compared to solutions 
obtained by the GUST3D frequency-domain solver both on the airfoil surface and in the 
far field. Grid density and domain size studies are included. 
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1. Introduction 


Computational aeroacoustics (CAA) is concerned with the time-accurate solution of 
flow and acoustic phenomena over long periods of time. In the problems of interest, the 
flow has a range of length and time scales, and bodies with complex geometry are in the 
computational domain. To accurately compute these unsteady problems, high-order 
finite-difference schemes and optimized schemes have been developed, as well as high- 
accuracy boundary conditions. For an overview of the CAA field, the reader is referred to 
two excellent review papers' ^ 

One problem of interest is the noise radiated when a vortical gust impinges on an 
airfoil with realistic geometry. This problem appears in helicopter noise as well as in 
turbomachinery noise. In the past, this problem has been solved numerically by finite- 
difference frequency-domain solvers by Scott and Atassi^ "*. In the time domain, it has 
been solved for flat plates using linearized potential by Hariharan et. aF, using the 
linearized Euler equations by Tam et. al.'^, and using full nonlinear Euler equations by 
Lockard and Morris’. The nonlinear Euler work of Lockard and Morris has also been 
extended to NACA 0006 and 0012 airfoil geometries, both lifting and nonlifting^ 

In past work, the gust response of two Joukowski airfoil geometries has been 
computed and compared to the GUST3D solutions using the NASA GRC PFC6 code.®’^° 
The ability of the PFC6 solver to accurately capture the pressure disturbance distribution 
on the airfoil surface has been demonstrated in that work. 

The present work is a more in-depth investigation of the performance of the PFC6 
code, comparing results from the near- and far-field as well as on the airfoil surface. The 
results show the ability of the code to accurately predict surface as well as far-field noise 



radiation. Grid density and domain studies show the effects of the boundary conditions 
and grid spacing on the solution obtained. 


2. Mathematical and Numerical Formulation 


In this work, the 2-D nonlinear Euler equations are solved. In Cartesian coordinates 
these equations are written as: 
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Since the Joukowski airfoil has a complex geometry that does not lend itself to 
Cartesian grids, the equations were recast in generalized curvilinear coordinates. From 
previous numerical tests”, the chain-rule formulation was chosen as the most accurate 
form of the equations. The chain-rule curvilinear Euler equations are written as: 
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The code used in this work was the NASA GRC PreFactored Compact 6 (PFC6) 
code, which is a predecessor to the current NASA GRC BASS code. The time stepping 
method used by PFC6 was the low storage fourth-order nonlinear extension of Hu’s 5-6 
Low Dispersion and Dissipation Runge-Kutta scheme'^ by Stanescu and Habashi”. 

The spatial derivatives are calculated using the prefactored sixth-order compact 
scheme and explicit boundary stencils of Hixon.'^' At block boundaries, an 11- point 
explicit stencil was used. A 10th order explicit filter'^ was used at every stage of the 
Runge-Kutta solver to provide dissipation. The code is written in Fortran 90 to take 
advantage of the improved memory management and data structures of Fortran 90 as 
compared to Fortran 77. This scheme and code has been extensively validated on 
benchmark problems, performing very well (e.g.. Refs. 9,10,14,16-18,21). 



The grid used employed a C-grid topology (Figure 1), extending at least five chord 
lengths away from the airfoil in each direction. Several tests were conducted for larger 
domains which extended 10 and 20 chord lengths away. The smallest grid was 433 x 125 
(54125 total) points, with the middle grid being 605 x 240 (145200 total) points, and the 
largest grid having 790 x 271 (214090 total) points. All grids were generated using the 
commercial package GridPro^®. The grid was clustered algebraically in the normal 
direction (An = 0.01) and near the trailing edge point (Ax = 0.01), as shown in Fig. 2. A 
stretching ratio of 1.05 was used to a far-field spacing of Ax = Ay = 0.106. This spacing 
gave a minimum of 15 points per gust wavelength. While this is far more than the six 
points per wavelength that the scheme could use on a uniform Cartesian grid, the dense 
grid was used to ensure that resolution problems would be minimized for this test. 

For the calculation, the grid was generated as a single block and the domain was 
distributed by the code for parallel computations on an arbitrary number of processors. 
For this 2D test case, tests on both a cluster of SGI Octane workstations as well as on a 
cluster of Apple Power Macintosh G4/450 machines showed excellent parallel 
performance for up to 15 machines. 


3. Problem Description 

In this set of problems, a simple-harmonic vortical gust convects past a 12% thick 
Joukowski airfoil. The gust has the distribution: 
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Here, e = 0.02. 

Two gust configurations were tested. The first was a transverse gust: 


a = 2k 


P = 0 (9) 

co = 2kM 


and the second was a two-dimensional gust: 


a = 2k 

P = 2k ( 10 ) 

co = 2kM 


where M is the freestream Mach number. 

For each gust distribution, two reduced frequencies were tested: k = 0.1, and k = 1.0. 
Figures 3 and 4 show instantaneous v- velocity distributions for the cambered airfoil cases 
at k = 1.0, while Figures 5 and 6 show the instantaneous pressure distributions for those 



cases. Note that the two-dimensional gust propagates at a 45° angle to the airfoil. Note 
also the clean boundary treatments and lack of reflections for both vortical and acoustic 
waves. 

The mean flow at infinity is defined as; 
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where M = 0.5 and y = 1.4. 

Two airfoil geometries are tested in this work. The first geometry is a symmetric 
airfoil at a zero degree angle of attack. The second geometry is a cambered airfoil 
(camber ratio = 0.02) at a two-degree angle of attack. 


4. Initial and Boundary Conditions 


For both cases the flow was initialized to the mean flow with the vortical gust 


superposed: 



( 12 ) 


p{x,y,0) = p 

u{x,y,0) = u+u^^,{x,y,0) 
v{x,y,0) = v + v^^,{x,y,0) 
p{x,y,0) = p 


At the wall, Hixon’s inviscid curvilinear wall boundary condition^® was used, 
modified to set the normal momentum to zero at the wall at each Runge-Kutta stage. 

At the inflow boundary, the acoustic radiation condition of Tam and Webb^‘ was used 
on the outgoing perturbations. For example, the outgoing u-velocity perturbation Ug^ was 
defined as; 


^BC ~ ^boundary ^ ^gust ( 13 ) 

At the outflow boundary, Tam and Webb’s radiation outflow condition^* was used 
with no correction for the outgoing vortical gust. 

At the trailing edge point on the airfoil, the C-grid topology defines the airfoil 
geometry incorrectly, causing numerical inaccuracy. While an 0-grid geometry would 
define the trailing edge geometry properly, the C- grid was chosen due to the excessive 
number of grid points and the associated small time steps that the O- grid would require 
to accurately resolve the sharp trailing edge. 

At the trailing edge, an upper and lower wall condition is calculated; this makes the 
trailing edge point double-valued. The wake is single-valued, however; the upper wake 
block is allowed to set the wake line. The discontinuity in the boundary condition on the 
surface line as it enters the wake causes a loss of accuracy near the trailing edge, which is 



especially apparent for the lifting, high-frequency results. To reduce the effect on the 
global solution, points were clustered near the trailing edge as shown in Figure 2. Figure 
7 shows the effect of the discontinuity on the pressure contours near the trailing edge of 
the cambered airfoil. The effect was much stronger on the cambered lifting airfoil than on 
the symmetric nonlifting airfoil. 


5. Results 

The PFC6 code was run until the lift coefficient settled to a simple harmonic state, 
corresponding to a nondimensional time of 210, requiring 55 hours using the largest grid 
on a 15 -processor Macintosh G4/450 cluster. Later optimization of the code reduced this 
time to 19 hours. As expected, the nonlifting airfoil case converged faster; the lifting 
airfoil case results were still changing very slightly from one period of gust oscillation to 
the next. A time step of CFL =1.5 was used for all calculations, giving approximately 
911 time steps per cycle of vorticity for the k = 1.0 case and 9110 time steps per cycle of 
vorticity for the k = 0. 1 case on the smallest grid. 

The mean pressures on the airfoil for the 2-D gust cases are shown in Figure 8, 
compared to FL036 results The effect of the trailing edge condition is apparent in both 
figures; however, the effect is localized near the trailing edge. In both figures, some 
oscillations are seen near the peak of the pressure curve; this is due to the high gradients 
in both the flow properties and the changes in grid spacing near the leading edge. 
However, the effect on the mean solution is minimal. It is seen that the lifting airfoil peak 


pressures are consistently overpredicted; studies found that enlarging the computational 
domain caused the predicted lift to decrease. 

Figures 9 and 10 show the RMS pressure disturbance distribution on the airfoil. 
Again, some oscillations can be seen near the peaks of the pressure disturbance, but these 
were found to have little effect on the solution. The trailing edge condition has some 
effect on the solution on the airfoil, particularly for the high- frequency gusts. However, 
as shown in Figures 9 and 10, the code is predicting the changing pressure distributions 
for all eight test cases very well, accurately capturing the effects of airfoil geometry, gust 
geometry, and gust frequency. 

Figures 1 1 and 12 show the acoustic intensity at a distance of one chord length away 
from the centerpoint of the airfoil. Here the comparison is very good for the low- 
frequency cases and still good for the high-frequency cases. Notice that the agreement is 
better for the nonlifting cases than for the loaded airfoils; however, all comparisons are 
very acceptable. 

Figures 13 and 14 show the acoustic intensity at a distance of four chord lengths away 
from the centerpoint of the airfoil. The GUST3D solution at this location was obtained 
using a Kirchhoff approach.^^ The Kirchhoff surface was a circle of radius two chord 
lengths, centered about the airfoil center. Again, the comparisons are very good for the 
low-frequency cases, and the high-frequency results are acceptable. 



6. Grid Density and Domain Size Studies 


Two solution convergence studies were performed in this work. In the first study, the 
effect of grid density on the solution was tested for the low-frequency, ID gust case on 
the cambered airfoil. A dense grid was generated which contained 677 x 233 points 
(157741 total points), but still using the small domain (5 chord lengths to the outer 
boundary). In this way, the effect of both clustering more grid points at the high-gradient 
leading edge region and increasing the number of grid points per wavelength could be 
determined. As Figure 15 illustrates, there was very little difference in the computed 
solutions for both grid densities. Thus, all other cases were run using grids closely 
resembling the small- domain grid near the airfoil. 

The next study focused on the effect of the boundary proximity to the airfoil. For the 
cambered, lifting airfoil case in particular, the distance to the far-field boundary condition 
had a large effect on the computed results. Figures 16-18 show the effect of increasing 
domain size on the acoustic intensity results at one chord length away from the airfoil, 
with the domain boundaries ranging from 5 to 20 chord lengths away from the airfoil. It 
was found during these numerical investigations that the outer boundary effect was 
minimal on the near-field solution as long as the outer boundary was located at least 10 
chord lengths away from the airfoil. 




7. Conclusions 


A test matrix of eight gust-airfoil interaction problems were solved using a 
parallelized prefactored sixth- order compact scheme with 10th order filtering, run on 
distributed computer clusters. These problems tested the accuracy of the code on 
stretched, curvilinear grids with nonlinear flows. In all cases, the code was robust and 
converged well. In these tests, the code predicted the trends and magnitudes of the airfoil 
gust response very well for the eight test cases, while showing some effect from the 
boundary condition reflections and the close proximity of the computational boundaries. 
The effect of larger domain sizes was investigated, and it was found that the near-field 
solution became invariant once the domain boundary was at least 10 chord lengths away 
from the airfoil. A grid density study was also performed, and it was found that the 
solution was effectively grid-independent for these calculations. 
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Figure 1: C-Grid Used for Joukowski Airfoil (433 x 125) 
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Figure 2: Closeup of Symmetric Airfoil Grid 







Figure 5: Instantaneous Pressure Distribution 
(k = 1.0, ID gust, cambered airfoil) 



Figure 6: Instantaneous Pressure Distribution 
(k = 1.0, 2D gust, cambered airfoil) 
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Figure 7: Effect of Trailing Edge Grid Singularity on Instantaneous Pressure 
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Figure 8: Mean Pressure Distribution on the Airfoil Surface 
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Figure 10; RMS Pressure Disturbance Distribution on the Airfoil Surface 
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Figure 11: Acoustic Intensity Distribution at R = 1.0 

(k = 0.1) 
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Figure 13: Acoustic Intensity Distribution at R = 4.0 
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